#=
=#


include("../src/qmcmary.jl")
using ..qmcmary

using Random
using Dates


"""保存hs场的位型"""
function savehs(hscfg, head, fname)
    fil = open(fname, "a")
    write(fil, head)
    write(fil, "\n")
    Nt = size(hscfg)[1]
    for ni in Base.OneTo(Nt)
        write(fil, string(hscfg[ni, :]))
        write(fil, "\n")
    end
    close(fil)
end



function forward(ehk, Nt, Nx, Ng, Ui, psi, the)
    nx = sin.(psi).*sin.(the)
    ny = cos.(psi).*sin.(the)
    nz = cos.(the)
    println(ny)
    println(nz)
    fname = "tmp"
    fil = open(fname, "w")
    write(fil, string(psi)*"\n")
    write(fil, string(the)*"\n")
    close(fil)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, ehk, hscfg, Ui, nx, ny, nz; Z2=true)
    for _ in Base.OneTo(100)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, ehk, hscfg, Ui, nx, ny, nz, sslen, bmats, allbmats, ss; Z2=true)
    end
    nup_i = zeros(2*Nx)
    tph = 0.0
    nsp = 100
    for idx in Base.OneTo(nsp)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, ehk, hscfg, Ui, nx, ny, nz, sslen, bmats, allbmats, ss; Z2=true)
        tph += ph
        savehs(hscfg, string(idx), fname)
    end
    println("forward sgn: ", tph/nsp)
end


function step_theta1(ehk, Nt, Nx, Ng, Ui, psi, the)
    fid = open("tmp", "r")
    psi2 = readline(fid)
    println(psi, psi2)
    the2 = readline(fid)
    println(the, the2)
    #
    nx = sin.(psi).*sin.(the)
    ny = cos.(psi).*sin.(the)
    nz = cos.(the)
    #println(ny)
    #println(nz)
    #
    #重新读取cfg
    ttbar = zeros(Nx)
    tpbar = zeros(Nx)
    tph = 0.0
    bistr = "0"
    println(Dates.now())
    while !eof(fid)
        bistr = readline(fid)
        #println(bistr)
        hscfg = Matrix{Int}(undef, Nt, Nx)
        for bi in Base.OneTo(Nt)
            cfg = readline(fid)
            cfg = split(cfg[2:end-1], ',')
            for xi in Base.OneTo(Nx)
                hscfg[bi, xi] = strip(cfg[xi]) == "1" ? +1 : -1
            end
        end
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, ehk, hscfg, Ui, nx, ny, nz; Z2=true)
        ###
        gf, ph = eq_green_scratch(ss)
        #println(ph)
        tph += ph
        #
        nxbar, nybar, nzbar = meas_grad(ss, allbmats, ehk, hscfg, Ui, nx, ny, nz; Z2=true)
        pnxpt = sin.(psi) .* cos.(the)
        pnypt = cos.(psi) .* cos.(the)
        pnzpt = -sin.(the)
        pnxpp = cos.(psi) .* sin.(the)
        pnypp = -sin.(psi) .* sin.(the)
        pnzpp = zeros(Nx)
        println("nzbar", nzbar)
        println("nybar", nybar)
        println("nxbar", nxbar)
        thetabar = @. nxbar * pnxpt + nybar * pnypt +  nzbar * pnzpt
        ttbar += real(thetabar)
        psibar = @. nxbar * pnxpp + nybar * pnypp
        tpbar += real(psibar)
        if bistr == "100"
            println(Dates.now())
        end
    end
    close(fid)
    nsp = parse(Int, bistr)
    println(tph/nsp)
    return ttbar/nsp, tpbar/nsp
end



function optimal(L)
    #
    Nx = 3*L^2
    psi = 0.25*π*ones(Nx)
    the = 0.25*π*ones(Nx)
    println(psi)
    println(the)
    hk = lattice_kagome(ComplexF64, L, -1.0+0.0im)
    Ui = 4*ones(Nx)
    Nt = 60
    Ng = 4
    for ts in Base.OneTo(100)
        #Random.seed!(678231)
        forward(hk, Nt, Nx, Ng, Ui, psi, the)
        tb, pb = step_theta1(hk, Nt, Nx, Ng, Ui, psi, the)
        the += 0.1*tb
        psi += 0.1*pb
        println("t", tb, pb)
        println(the, psi)
    end
end


optimal(9)

